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Abstract. The acceleration of ultrahigh energy nuclei in fast spinning newborn pulsars can 
explain the observed spectrum of ultrahigh energy cosmic rays and the trend towards heavier 
nuclei for energies above 10^^ eV as reported by the Auger Observatory. Pulsar acceleration 
implies a hard injection spectrum (~ E~^) due to pulsar spin down and a maximum energy 
ErasLx ~ Z 10^^ eV due to the limit on the spin rate of neutron stars. We have previously shown 
that the escape through the young supernova remnant softens the spectrum, decreases slightly 
the maximum energy, and generates secondary nuclei. Here we show that the distribution of 
pulsar birth periods and the effect of propagation in the interstellar and intergalactic media 
modifies the combined spectrum of all pulsars. By assuming a normal distribution of pulsar 
birth periods centered at 300 ms, we show that the contribution of extragalactic pulsar births 
to the ultrahigh energy cosmic ray spectrum naturally gives rise to a contribution to very 
high energy cosmic rays (VHECRs, between lO^*^ and 10^* eV) by Galactic pulsar births. 
The required injected composition to fit the observed spectrum depends on the absolute 
energy scale, which is uncertain, differing between Auger Observatory and Telescope Array. 
The contribution of Galactic pulsar births can also bridge the gap between predictions for 
cosmic ray acceleration in supernova remnants and the observed spectrum just below the 
ankle, depending on the composition of the cosmic rays that escape the supernova remnant 
and the diffusion behavior of VHECRs in the Galaxy. 
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1 Introduction 

The origin of cosmic rays continues to challenge our understanding after a century of ob- 
servations. Space and balloon based observatories have made precise measurements of the 
spectrum and composition of cosmic rays up to < 10^^ eV per particle. Above these energies, 
observatories on the ground have studied extensive air showers up to 10^*^ eV. The bulk of 
the cosmic ray flux is believed to be accelerated in Galactic supernova remnants (SNR) [1-3] . 
This long held notion fits well the observed spectrum up to 10^^ eV [4]. Above these energies 
a new component is needed to explain the spectrum and observed composition. This new 
component may be Galactic, as suggested in [5, 6], or extragalactic as proposed in [7, 8]. 
Constraints derived from the lack of strong anisotropics in cosmic ray arrival directions and 
the observed structure of the Galactic Magnetic field show that above 2 x 10^^ eV cosmic rays 
are extragalactic [9]. However, the transition from Galactic to extragalactic is expected to 
occur well below 10^^ eV, with models spanning the very high energy (VHE) range between 
10^^ eV and 10^^ eV with "dip" models around 10^'' eV [7, 8] and "ankle" transition models 
around 10-"^^ eV (see, e.g., [10]). 

The study of ultrahigh energy cosmic rays (UHECRs), from 10^^ eV to 10^'' eV, has pro- 
gressed significantly with the advent of giant airshower observatories such as the 3,000 km^ 
Pierre Auger Observatory in Mendoza, Argentina [11] and the 700 km^ Telescope Array (TA) 
in Utah, USA [12, 13]. The spectrum, sky distribution of arrival directions, and composi- 
tion indicators are well measured over a large range of energies. Difl'erences in reports from 
the two major observatories include a 20% shift in absolute energy scale (-EAugcr — O.SE'ta) 
and the diff'ering trends of composition indicators at higher energies. Currently the most 
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extensive dataset on composition indicators, such as the average and the RMS of the depth 
of shower maximum (Xniax)j has been pubhshed by the Auger collaboration and shows a 
departure from a composition consistent with lighter nuclei at 10^*^ eV to a trend towards 
heavier nuclei above 10^^ eV [14, 15] (see Fig. 2). TA reports shower behaviors consistent with 
protons [16]. The discrepancies in composition reports and the difference in absolute energy 
scale make it difficult to constrain proposed models for the origin of UHECRs. Fortunately, 
a cross-experiment effort to understand these discrepancies is currently on-going. 

Here we show that the fast spinning pulsar birth model described in [17, 18] can explain the 
observed spectrum (both the Auger and the TA spectra) and the composition trend observed 
by Auger [14, 15]. To fit these two observables we allow the freedom to vary the percentage 
of different elements that are accelerated in the pulsar wind divided into 3 groups: Hydrogen, 
Carbon group (CNO), and Iron. Although the surface of the rotating neutron star is thought 
to be a natural source of Iron, X-ray spectra of pulsars have shown evidence for Helium 
[19], and Carbon, Oxygen, and Neon [20, 21]. In addition, the material leftover from the 
progenitor star is likely to be in the Carbon group (mostly Oxygen, see e.g., [22-25]). As we 
show below, to fit the Auger spectrum a balanced ratio between Hydrogen and CNO, with 
a less fraction of Iron suffices, while to fit the TA spectrum a higher percentage of Iron is 
needed since the spectrum extends to higher energies (see Fig. 1). The composition selected 
by the Auger spectrum gives a good fit to the Auger average shower maximum ((Xmax)) and 
the fluctuations around the mean (RMS(Xinax)), see Fig. 2. 

The contribution of young pulsars to the Galactic cosmic ray flux is related to the extra- 
galactic contribution through the distribution of birth parameters. Here we calculate this 
contribution for different UHECR scenarios. The end of the Galactic spectrum is highly 
dependent on Galactic diffusion parameters and the history of the most recent pulsar births 
in the Galaxy (such that the flux in the VHE region varies over long time scales). However, 
the differences in the UHECR data, such as the different absolute energy scales of the Auger 
and TA observatories, give rise to significant changes in the predictions for the Galactic com- 
ponent. The Auger fit implies Galactic pulsars as the main contributors to VHECRs from 
10""^^ eV, while with the TA fit, they start dominating from 10^ ''■^ eV (see Fig. 4). 

As an acceleration model that naturally produces heavier nuclei at ultrahigh energies, the 
fast pulsar birth model illustrates the tension in fitting the UHECR spectrum, composition, 
and anisotropy measurements over the highest two decades in energy (for further discussion 
on this tension see, e.g. [10, 26-32]). The relevant energy range is bracketed between a 
transition from Galactic to extragalactic and the cosmological GZK propagation effect. The 
observational uncertainties are still too large to allow for a clear picture of the source re- 
quirements. Great progress will occur when the leading experiments at the highest energies, 
Auger and TA, reach a better agreement on the absolute energy scale and the composition 
as a function of energy. In addition, increase statistics in observations, e.g., from the JEM- 
EUSO space mission [33], can reveal a clearer anisotropy picture, which will greatly improve 
our ability to zero in on the best model for the origin of ultrahigh energy cosmic rays. 
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2 Newborn Pulsars as cosmic ray sources 



The acceleration of particles in pulsar environments has been suggested since their discovery 
[34]. While nearby pulsars show direct evidence of accelerated electrons and positrons, the 
acceleration of hadrons is still unclear. The suggestion of pulsars as cosmic ray accelerators 
of VHECRs has been discussed in [35-39] while for UHECRs the main proposals are by 
[17, 40]. In [17] Iron nuclei stripped off the neutron star surface are accelerated to UHEs 
by the fastest spinning young neutron stars with typical pulsar magnetic fields (between 
10^2 and 10^^ G). They derived the maximum energy and spectrum due to the spin down 
of young pulsars (J oc E~^) that we use below. Neutron stars with much larger surface 
magnetic fields, i.e., magnetars, have been proposed as sources of ultrahigh energy protons 
by [40]. Given their faster spin down rate, the acceleration to UHEs is at earlier stages of the 
pulsar evolution when gravitational radiation is significant and a disruption of the supernova 
envelope is needed to allow the escape of accelerated particles [40, 41]. 

In [18] we discussed the acceleration and escape of UHE cosmic rays from newly born pulsars. 
Here we follow the same notation describing a pulsar with magnetic dipole moment fi = 
A*30.5 lO^'''^ cgs and rotational speed ^1 = 0.4 10^ s~^, that accelerate particles of charge Z = 
26^26 up to the energy [17]: 

E{f,, n,z)^9x 10^0 m m.5 ev (2.1) 

where < < 1 is the pulsar acceleration efficiency, that is, the fraction of the full potential 
drop a charged particle experienced during its way out of the pulsar wind. We take 7^3 = 77/0. 3 
in this work. 

As the pulsar spins down, the energy of the particles produced in the wind decreases. A 
particle of energy E is accelerated in the wind at time: 



9 X 10^0 eV\ ^26773^5 



tspi„(^) ~ 3 X loM ) ""''"^"^ s , (2.2) 

V / A^30.5 

neglecting energy losses via gravitational waves (see [40, 41]). 

The cosmic ray spectrum injected over time by a pulsar with inertia / = 145 10^^ gcm^ can 
then be calculated, assuming the Goldreich- Julian charge density, ncj = B i^/iirZec, [42] 
is entirely tapped in the outflow [17]. The total charge density is the sum of the densities 
for each chemical group uq] = fz nQ]{Z), where fz is the fraction of the Goldreich- 
Julian charge density in particles with charge Z injected into the pulsar wind and iz satisfies 
^^f^ = 1. The injected spectrum of a given species Z is then given by: 

^ (/i, n,Z) = 5x 10^3 iz 745(^26 /U30.5 ^2o)"' eV"^ (2.3) 

In the pulsar model, the more energetic particles are accelerated at earlier times, when the 
surrounding supernova envelope is denser making the escape more difficult. Assuming a 
supernova expands with constant speed Vcj and explosion energy E'cj = -E'cj,52 10^^ erg, the 
mean density of the supernova ejecta with mass Mcj^io = M^j 10 is: 



-3- 



Numerical results of the escaped spectrum can be found in Figure 3 of [18] for Hydrogen 
and Iron. Here we also include the results of similar calculations performed for the escape 
of injected Helium and CNO. The traversal of the supernova ejecta exponentially cuts the 
injection spectrum at: 



i^eut,Z ^ 2.3 X 1019 Zir?3/45/^3"o'5^4T,\o^Jj(52 (^) (2.5) 
3.6 X 1020Z26r?3/45/^3"a5^e7io^S ( v) ' ^^'^^ 

where a is the cross section of the hadronic interaction between cosmic rays and ejecta par- 
ticles. In the relevant interaction energies, for proton-proton interactions it is approximately 
(7p ~ 130 mb while it reaches apc ~ 1-25 b for Iron-proton interactions (note that we used cross 
sections with full energy and composition dependence in the simulations). For E < Ecut,z, 
the spectrum is not affected by the interactions. Above -Ecut.Zi the escaped spectrum can be 
approximated by: 

^ (/., ^,Z) = ^ (M, n, Z) e^-^/^- (2.7) 

The chemical composition of the ejecta after steady and explosive burning of supernova II-P, 
lb, and Ic has been studied by a number of authors, e.g., [22-25]. Core-collapse supernova 
progenitors have large abundances of Oxygen produced by Helium burning during the life of 
the massive star, and part of this Oxygen is burnt during the explosion, which produces Si, 
So, Mg, Ca and Ni. Note that in most cases, Hydrogen is the dominant component in the 
envelope to consider for the escape of cosmic rays, as studied in [18]. The effect of escape from 
SN envelopes which are not proton dominated is small once the summing over all pulsars is 
considered as below. 

The source of UHECRs in our model are the rare, extremely fast spinning, young pulsars. 
The majority of pulsars were born spinning slower and therefore contribute to the flux of 
lower energy cosmic rays. According to [43], the distribution of pulsar birth spin periods, 
f[P = 27r/il), is normal, centered at 300 ms, with standard deviation of 150 ms. We 
assume that proto-pulsars spinning initially below the minimum spin period of a neutron 
star Pmin — 0.6 ms [44] evolve to a stable figuration with P = Pmm- Hence we assume 
/(-Fmin) = Ylpf{P < -fmin)- The initial magnetic field follows a log-normal distribution 
/(^) with (log(B/G)) - 12.65 and cJiogB ~ 0.55 [43]. 

The total cosmic ray spectrum contributed by the entire pulsar population is then 

^ (^) = / ^ if^, Z) /(/.) f{^) df, dn (2.8) 

Cosmic rays injected by Galactic and extragalactic (EG) pulsars travel through the Inter- 
stellar Medium (ISM) and the Intergalactic Medium (IGM) on their way to Earth. The 
corresponding propagation effects also affect the observed spectrum and composition and are 
discussed next in Section 3 and Section 4. 
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3 Propagation from Extragalactic sources 



On average, the pulsar birth rate in the Galaxy is z/g — 1/60 yr [45]. For a galaxy density 
of ngai — 0.02 Mpc~^ and an energy loss time for cosmic rays with energy above 10^^ eV of 
cTioss ~ 792 Mpc [46], the average flux of cosmic rays from EG sources can be estimated to 
be: 



diVEG 



dEdtdAdn 



dN I 

dEA^""'"'' 
dN , 

TT^ /s 7 X 10" 
dE ■' 



eV 



m 



sr 



(3.1) 
(3.2) 



where /s is an overall factor used to fit the model prediction to the measured UHECR flux. 
In the pulsar model, it can be interpreted as the fraction of total flux of pulsar births required 
to account for the observed flux of UHECRs out of the total pulsar birth rate, /s < 1 can 
be due to efficiency factors such as variations in the core-collapse geometry, poorer injection 
efficiency, or a lower hadronic density in the pulsar wind than the Goldreich-Julian density. 

In order to estimate the observed spectrum of UHECRs, after their propagation through 
the extragalactic medium, we rescaled the simulation output of [47] from 10^^ eV up to 
the maximum acceleration energy, -Emax by the injected spectrum of pulsars. The effect of 
the well-known Greisen-Zatsein-Kuzmin interactions of UHECRs with the cosmic microwave 
background [48, 49] (as well as other cosmic photon backgrounds) is included in these cal- 
culations, however, in the pulsar model, the GZK effect is secondary to the effect of -Emax 
which is set by the fastest spin rate that neutron stars can reach J^max lO^'-^s"^ [44]. 

The UHECR spectrum and composition ratios for a given pulsar, after escape through the 
surrounding supernova ejecta, were calculated in [18]. At the highest energies, the spectrum 
is mostly determined by the superposition of exponential cutoffs up to -Emax; coming from 
the effect of escaping through the supernova ejecta. In order to account for the whole range 
of pulsar spins and magnetic fields, we ran 19 x 19 simulations of sets of (P, log/i), and 
integrated the obtained spectra over the overall pulsar population following Eq. 2.8. 

Figure 1 shows the propagated energy spectrum of cosmic rays from EG pulsars. In the 
middle and bottom panels, the pulsar emissivity is assumed to follow the star formation rate 
computed by [50]. Note that a star formation history following the model of [51] only changes 
slightly the results. For comparison, we also present a case where the source emissivity 
remains constant over time in the top panel, which we call Auger-uniform. (For a discussion 
on the influence of the source evolution model on the UHECR spectrum, see [47].) In the 
Auger-uniform case, the sources were assumed to be uniformly distributed in space, normally 
distributed in log /Li and P as in Equation 2.8 and with pulsar wind acceleration efficiency 
r] = 0.3. 

The overall normalization factor /s 0.05 for the Auger-uniform (top panel) and Auger 
cases (middle panel), and /s ~ 0.08 for the TA case (bottom panel), to take into account 
the difference in energy scale between the two experiments. In the Auger-uniform case the 
injection composition is 50% H (fn = 0.5), 30% CNO (fcNO = 0.3) and 20% Fe (fpe = 0.2). 
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In the Auger case fn = 0.65, fcNO = 0.20, and fpc = 0.15. These composition ratios were 
chosen to fit the Auger spectrum [52] and composition [14, 15]. In the TA case, an injection 
of fu = 0.5 and fpc = 0.5 provides a better fit to the Telescope Array (TA) spectrum [53]. 

In the Auger case, our model over-produces the total flux between lO^*^'^ — 10^^ eV. In 
this study we do not include the possible effect of a magnetic horizon as in [8, 54-56] (see 
Section 6). This effect would harden the extragalactic spectrum around the ankle region, and 
thus enable a better spectral fit with the star formation rate scenario. It is also possible that 
the extragalactic pulsar population that produces UHECRs has a biased evolution (see [18]) 
that can be weaker compared to the classical star formation rate. One would then obtain 
a better fit to the spectrum, as in the Auger- uniform case, but, as discussed below, to the 
expense of a good fit to the composition around the ankle region. 

The mean atmospheric depth (Xmax) and its fluctuations, RMS((Xmax)), of the Auger case 
compared with Auger measurements [14] are presented in Fig. 2. Four hadronic interaction 
models EPOSvl.99 [57], QGSJETOl [58], QGSJETII [59] and SIBYLL2.1 [60] were used in 
this calculation, which gives the range of the blue shaded region (also see [61] for fits to the 
Auger data with different hadronic interactions). Our estimates include the contribution of 
Galactic pulsars, as calculated in the following section. 

Our results fit the observations better with the EPOS model. (Note that we used the EPOS 
model in our simulation of the escape of accelerated particles from the supernova ejecta.) 
The Auger and Auger-uniform cases follow well the trend of a transition from light to heavy 
elements as measured by the Auger Observatory. In the Auger case, the composition is 
dominated by the contribution of extragalactic pulsars down to ~ 10^^ eV. In the Auger- 
uniform case however, pulsars become an underdominant source below E ~ 10^^-^ eV, as the 
extragalactic contribution cuts off at the ankle, and the Galactic pulsar population provides 
only a fraction of the total flux below. As a result, the composition below the ankle appears 
heavier than measured, and another light component is needed (which raises some questions 
about the origin of this other component as discussed at the end of the Section 5). Under our 
chosen assumptions, the TA case does not fit the constant light composition measurements 
of TA and is a bit heavier than the Auger composition. 

Because protons dominate the flux at ankle energies, and Helium nuclei are mostly dissociated 
into protons during the propagation, the injected protons can be mostly interchanged to 
Helium without affecting the spectrum signiflcantly. For example in the Auger-uniform case, 
the composition and spectrum after propagation remains nearly unchanged for an injection 
of fn = 0, fHc = 0.5, fcNO = 0.3, and fpc = 0.2 (Figure 3). 

4 Propagation from Galactic sources 

For the propagation of cosmic rays accelerated by Galactic pulsars, we model the turbulent 
Galactic magnetic fleld as a cylindrical halo of radius i?Gai = 15kpc, of height above (or 
below) the Galactic plane typically ~ 2 — 8kpc [62], of coherence length Ic = 10 — 100 pc 
and strength B = 3 fiG (see [63] and references therein) . The Larmor radius of a particle 
reads tl = 13.8 Eis Z:^^ {B/3 ^lGy^ pc. 
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The flux of cosmic rays accelerated by Galactic pulsars can then be calculated as 



dEdtdAd^ dE 47r Vbai V ^ ' 

where /s is the same scale factor introduced for the EG component, is the average birth 
rate of pulsars in our Galaxy, and Vcai = ^vr HRq^^ is the volume of the Milky Way. The 
average number of sources in the Galaxy that contributes to the observed spectrum at a 
given energy E is Ns{E) = (5iobs(-S)/5tbirth (-£■), where 6tobs{E) is the timescale over which a 
source can can contribute to the observable cosmic rays of energy E, and (^tbirth(-E')~^ is the 
birth rate of Galactic sources that can produce particles of energy E. Following a Poisson 
distribution, the probability that currently at least one source is contributing in the Galaxy 

ig ^ _ g — '5*obs/'5ibirth ^ 

If the scattering length of the particle (distance over which its deflection angle becomes 
69 ^ 1) is shorter than the height of the Galaxy, the propagation will be mostly diffusive, 
and the source observation time is equivalent to the particle escape time from the Galaxy: 
^iohs = ''"esc- At energies above the knee (~ 10^^ eV), which is the main concern of this paper, 
nuclei spallation is negligible [64] and the nuclei escape time can simply be estimated with 
the Leaky box model 

T,sc{E,Z,k) = —, (4.2) 

where D is the diffusion coefficient in the Galactic magnetic field, that can be estimated 
empirically as in Equation A2 of [65] 



D{E,Z,l,)^Dori^c 



. '2/3' 



(4.3) 



The coefficient a depends on the turbulence and structure of the Galactic magnetic field. The 
normalization Dq is set at energies where particles are in the Kolmogorov regime (rL <S /c)) 
using measurements of the boron to carbon ratio in our Galaxy. We follow the estimates of 
[64]: 

D{R) = 1.33 X 10^8 //kpe f ^ j cm^s-^ , (4.4) 

where -ffkpc = i?/lkpc, and the particle rigidity R = E/Z. Notice that in the diffusive 
regime, the diffusion coefficient scales as E~^ , with (5 = 1/3, a value which seems to be 
favored by observations, e.g., [64]. A larger value of /3 would imply a faster escape out of the 
Galaxy of particles at very high energies. For example, (3 = 0.6 would imply that all particles 
above 4 x 10^^ GV travel rectilinearly, which is inconsistent with anisotropy measurements 
[66-68]. 

For scattering lengths larger than iJ, the propagation is quasi-rectilinear. In principle, the 
observation time 5tohs is then equivalent to the dispersion of particle arrival times at = 
DgdO"^ /{4:c), where Ds is a typical source distance [69-71]. The transition from a totally 
diffusive regime to a quasi-rectilinear regime can be modeled using Eq. 4.3 with a = 1. 
For reasonable assumptions on the structure of the magnetic field, one can then assume 
Stohs = o-t ~ Tcsc at high rigidities. 
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5 Very High Energy and Ultrahigh Energy cosmic rays from Pulsars 



Assuming that cosmic rays injected in the ISM by Galactic pulsars have the same composition 
as those from EG pulsars, we can calculate the contribution of both Galactic and extragalactic 
pulsars. Figure 4 show the results for the Auger, Auger-uniform (top and middle panels) 
and the TA cases (bottom panel) where we also show the total energy spectrum and the 
spectrum decomposed into three components (Hydrogen, intermediate, and Iron) as reported 
by KASCADE [74], KASCADE-Grande [72] and the Auger and TA spectra. The EG spectra 
(in dash lines) in Fig. 4 are the same as those in Fig. 1. 

The Galactic spectrum of each element group (in solid lines) has three regimes. For exam- 
ple, if we consider the Iron branch in the Auger-uniform case (top panel), it shows a totally 
diffusive behavior between E = 10^^'^ eV and 10^'^'^eV (where tl = Ic)- In this range the 
Galactic propagation softens the intrinsic spectrum by 1/D{E) ~ E~^'^ Z^'^. The second 
regime lies roughly between 10^^'^ eV and 10^^ eV, where tl > Ic and particles random walk 
with small deflections leaving the Galaxy a bit faster, so D{E) ~ E'^ Z~'^ (see Equation 4.3). 
However the second regime is overwhelmed by the last component, which comes in above 
10^^ eV when particles escape faster than the rate they are born, resulting in an event prob- 
ability of (1 — e~''*°'^=/^*'''"'*'^) and a time dependent flux. The three branches with Z = 1,7, 26 
for Hydrogen, CNO, and Iron have similar behaviors along the energy axis scaled by Z, 
i.e., same behavior for the same rigidity. The relationship between the amplitudes of the 
flux for each component combines the injection and propagation. The injected spectrum 
is inversely proportional to the charge dN^sc/dE oc Z~^. If each of the three elements are 
injected with fractions F{Z), then dN/dE oc Z~^''^ F{Z) in the totally diffusive regime and 
dN/dE oc Z^ F{Z) in the small deflection regime. 

Figure 5 contrasts the composition predictions of our models with the measurements of In A 
from Tunka ([76], filled square), Yakutsk ([77, 78], open square), C ASA-BLANC A ([79], 
open left triangle), HiRes/MIA ([80], open circle), HiRes ([81], filled downward triangle), 
KASKADE-GRANDE([72], filled diamond), Auger ([15], filled circle) and TA([73], open up- 
ward triangle) based on hadronic interaction model EPOSvl.99 [57]. The composition trends 
of our model reproduces some of the Galactic-Extragalactic transition features at the ankle 
as shown by the Auger and the KASKADE measurements. The contribution below the ankle 
depends on the dominant component in this energy range (defined here as more than 80% 
flux contribution) which may or may not be the pulsars (see e.g the Auger- uniform case). If 
the additional flux comes form acceleration in supernova remnants the composition is likely 
to be mostly Iron below the ankle (e.g., [82]). 

The expected anisotropy signal at the highest energies (above Eqzk = 6 x 10^^ eV) was 
already discussed in [18]. Because extragalactic pulsars are difficult to detect beyond our 
Local Group, and because of their transient nature, no direct correlation should be found 
between the arrival directions of the most powerful events and active sources. The distribution 
of arrival directions in the sky should trace the large scale galaxy distribution with a possible 
bias [71, 83], if the deflections experienced by the particles are small. At energies well 
below Eqzk, the arrival directions of particles emitted in extragalactic sources should be 
isotropic whatever their composition, as the deflections experienced by protons are already 
large enough to prevent any clustering. 
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Nevertheless, Galactic sources producing protons above energies ~ 10 eV could produce 
substantial anisotropics. In the regime of diffusive propagation, the anisotropy signal in a 
given direction can be defined as 6^ = (V^ncR/ncR) {3D/c), where ncR is the cosmic ray 
number density in the Galaxy as measured at the position of the Earth. Reference [84] 
showed that by assuming a homogeneous distribution of sources in the Galactic disc, the 
small-scale anisotropy signal can be simplified to (see also [85, 86]) 



Assuming diffusion parameters H = 2kpc and Ic = 20 pc, a strong anisotropy above rigid- 
ity i? ~ 4 x 10^ ''V would be expected from the non-diffusive regime when particles travel 
semi-rectilinearly. Measurements around 10^* eV indicate that cosmic rays at these energies 
are mostly light (see, e.g., KASCADE-Grande data), and are distributed isotropically in the 
sky [66-68]. The combination of these measurements require that particles around 10^^ eV 
(below the ankle) be protons accelerated in extragalactic sources. Galactic sources injecting 
protons at these energies would indeed induce strong anisotropics. This is a difficulty for 
most "ankle" -transition models, for which the extragalactic component becomes quickly un- 
derdominant below -Eankie = 10^^'^ eV. This is the case in particular for the Auger-uniform 
model, in which the transition at i^ankic implies that another light extragalactic component 
is needed to bridge the gap between a Galactic (e.g., SNR) component and the extragalactic 
pulsar contribution. 

In the Auger and TA cases however, for energies roughly below 10^^ eV, Galactic pulsars 
have a negligible contribution to the anisotropy signal, as other Galactic sources (e.g. SNR) 
dominate; for energies above 10^® eV, extragalactic sources dominate, which should produce 
a globally isotropic signal. For the energy range in between. Fig. 4 shows that the cosmic 
ray flux is mostly composed of CNO and Iron nuclei, that produce less anisotropy as they 
have smaller rigidity than protons at the same energy. In addition, [84, 87] show that effects 
of stochasticity in the spatial and temporal distribution of sources and turbulent magnetic 
field could cause significant fluctuation to the intrinsic anisotropy. In general our results are 
consistent with the anisotropy measurements [66-68]. 

This study reveals the tension between the spectrum, the composition, the anisotropy signal 
of cosmic rays above 10^^ eV and the need to have more precise measurements of these 
observables. The Auger-uniform model, though it can fit successfully the observed spectrum 
and the composition, requires a second extragalactic component of light elements around 
10^^ eV. The Auger case fits successfully the composition and alleviates the anisotropy issue 
at 10^^ eV, but one has to invoke a strong magnetic horizon effect to harden the spectrum at 
-E'anklc- Finally, the TA case is satisfactory for spectrum and anisotropy, but does not seem 
to account well for the observed composition. 

6 Discussion 

The results that we provided here for the Galactic component, depend on our choice of pa- 
rameters. Our fits to the observables of low energy cosmic rays were performed by adjusting 
key parameters of the Galactic magnetic fields: the coherence length of its turbulent com- 
ponent, Ic, and the height of its halo, H. Our knowledge on these parameters remain poor 
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([88] for a review), and leaves room to other shapes and normahzations for the Galactic 
component calculated here. For example, for H higher (lower) than 2 kpc, the field would 
confine -ff/2kpc times more (less) nuclei, as dNc^i/dE oc H^/D{E) oc H. The Galactic 
flux would hence move up (down) by H/2 kpc. A larger (smaller) coherence length of the 
field would enable particles with higher (lower) energy to be confined, and would shift the 
Galactic component to higher (lower) energies. Note that this effect is non linear, according 
to Equation 4.3. Further uncertainties in the modeling of particle diffusion, as discussed 
extensively in Ref. [64], could also noticeably modify the shape of the Galactic spectrum. We 
also assumed a uniform density of sources in our Galaxy to estimate the average properties, 
but at the very high energies, cosmic ray observables depend on the particular history of 
pulsar births in the Galaxy. 

On the other hand, our results on the UHE component due to the extragalactic pulsar 
population can be regarded as fairly robust. For UHECRs indeed, the shape of the propagated 
spectrum depends mostly on the injected spectrum and composition, the source emissivity 
evolution (see e.g., [47, 89]), and possibly on the presence of nearby sources at the highest 
energy end. 

The effects of the extragalactic magnetic fields are not taken into account in the propagation 
of UHECRs. At energies around the ankle, the presence of a mean field of order 1 nG could 
make the trajectories of particles so diffusive that they would not be able to propagate over 
the Hubble distance. This "magnetic horizon" effect could harden the spectrum of ^ £'0.i-o.3 
around the ankle region (e.g., [8, 54-56]). This effect would flatten the bump in the ankle 
region that appears in Fig. 1 and Fig. 4, and lead to a satisfactory fit of our model to the 
data for a broader range of Galactic magnetic field parameters. 

Note that the normalization factor of fg ~ 0.05 that we require for a good fit, is a physically 
reasonable value. This factor accounts for pulsar wind injection efficiencies and for the 
fraction of the total pulsar population producing cosmic rays. A value of fg ~ 0.05 can mean 
that only 5% of the pulsars are required to have the right configuration to produce cosmic 
rays. Assuming the gaussian distribution of pulsar parameters proposed by [43], less than 
0.3% of these 5% are pulsars born with periods less than 6 ms, that can accelerate particles 
to above 10^^ eV (these estimates can be modulated by the injection efficiencies). 

The present results also depend on the distribution function of pulsar parameters. We nor- 
malized here the overall Galactic and extragalactic spectra based on the data at the highest 
energies, i.e., with the fastest-spinning pulsars corresponding to the tail of the distribution. 
Slight variations in the bulk of the population could hence impact the shape of the spectrum, 
and thus the composition. For example, if f{P) cuts at Pmin instead of having those spin 
faster piled up to Pmim to meet the observed flux level at UHE, the Galactic component in 
the Auger case would over-produce VHECRs measured by Kascade. 

Recently, Ref. [90] discussed that muons produced by hadronic and photopion interactions 
during the acceleration process could be accelerated significantly before their decay, producing 
an enhanced neutrino flux at high energies. The authors consequently claim that most 
UHECR candidate sources involving unipolar induction should be ruled out, as they would 
overproduce neutrinos compared to the IceCube limits, via this mechanism. However, UHE 
ions injected at low latitude in rotation powered pulsars or magnetars barely meet background 
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hadrons while they 'surf-ride' in the force- free magnetosphere [40, 91]. Furthermore, we 
pointed out in Ref. [18], that the radiation fields in the pulsar wind are too low to impact the 
acceleration of UHECRs (see also [92], [37]). Note also that Ref. [90] assumes an injection in 
i<;~2.3 |-jy ijjjgar acceleration in their calculations, while our induction model implies a harder 
injection in (Eq. 2.3). This harder injection should result in a lower neutrino flux at 
the highest energies, than estimated by Ref. [90]. It thus appears that the muon acceleration 
mechanism should be negligible in our framework. 

Our scenario, as for any viable scenario, assumes that no recent events in our Galaxy that can 
reach E > 10^*'^ eV, in order not to overshoot the observed spectrum, and not to produce a 
striking anisotropy pattern in the sky. If a nearby pulsar were born with high enough spin to 
reach > 10^^-^ eV, then we should see a surge in the cosmic ray flux after a delay of thousands 
of years (for a similar treatment of GRBs see [93, 94]). 

7 Conclusion 

We showed how fast spinning pulsars can explain the observed spectrum of UHECRs (both 
Auger and TA) and the composition trend described by the Auger collaboration. To fit 
these two observables a mixed composition of Hydrogen, CNO, and Iron needs to escape 
the young supernova remnants accelerated via the fast spinning pulsar winds. To fit the 
Auger spectrum a balanced ratio between Hydrogen and CNO, with a minor presence of Iron 
suffices, while to fit the TA spectrum a higher percentage of Iron is needed. Determining the 
absolute energy scale is thus an important goal for current observatories, as it would help 
select between possible explanations for the origin of UHECRs. 

The composition mixture chosen to fit the Auger spectrum, gives a very good fit to the average 
shower maximum ((A^max)) and the fluctuations around the mean (RMS(Xmax)) observed 
by Auger. This is a unique aspect of the pulsar model as most models of astrophysical 
accelerators of UHECRs do not explain these shower maximum data. This challenge has led 
to the notion that the change in (Xmax) and RMS(Xmax) is due to new physics in hadronic 
interactions at these energies, which are well above those reached by the Large Hadron 
Collider. 

Another aspect of this model worth highlighting is the connection between parameters needed 
to fit the extragalactic component and the presence of a Galactic component from Galactic 
pulsar births in the very high energy range (between 10^^ and 10^^ eV). In the estimates 
presented here, the Auger-uniform fit implies an under dominant contribution to the flux 
of VHECRs, while the Auger and TA fits suggest that Galactic pulsars could be the main 
contributors to VHECRs. 
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Figure 1. Propagated energy spectrum of UHECRs from newly born pulsar population with log /i and 
P normally distributed, and wind acceleration efficiency 77 = 0.3. Simulation results were normalized 
at 10^^ eV with /g ~ 0.05 for the Auger and Auger-uniform cases, and fs — 0.08 for the TA case 
(due to difference in energy scale). The spectrum of each group of propagated nuclei are shown as 
in the legend box. Top (Auger-uniform case): a mixed composition of 50% Proton (fn — 0.5), 30% 
CNO (fcNO = 0.3) and 20% Fc (fpc = 0.2) was injected to fit the Auger spectrum [52]. The source 
emissivity is assumed to be constant over time. Middle (Auger case): 65% Proton, 20% CNO and 
15% Fe nuclei was injected also to fit the Auger spectrum [52], but the source emissivity is assumed 
to be follow the star formation rate computed by [50]. Bottom (TA case): 50% Proton, 0% CNO and 
50% Fe nuclei was injected to fit the TA spectrum [53]. The source emissivity is assumed to be follow 
the star formation rate computed by [50]. 
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Figure 2. (X,nax) (left column) and RMS(Xniax) (right column) of the Auger data [14, 15] (black 
crosses) and simulation results of the Auger- uniform (top), Auger (middle) and TA (bottom) cases as 
in Fig 1 (blue shaded region is for where pulsars contribute more than 80% to the total flux, hashed 
region when they contribute less). Flux from both Galactic and extragalactic pulsars arc taken into 
account. Four interaction models, EPOSvl.99 [57], QGSJETOl [58], QGSJETII [59] and SIBYLL2.1 
[60] were used to estimate the range of (Xmax) and RMS(Xmax) as listed in the legend box. The red 
and dark blue lines correspond to 100% P and 100% Fe. 
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Figure 3. Propagated energy spectrum (top), (Xmax) (middle) and RMS(Xniax) (bottom) of an 
alternative Auger-uniform case, with Proton interchanged to Helium at injection. A mixed compo- 
sition of 50% Helium (fne = 0.5), 30% CNO (fcNO = 0.3) and 20% Fe (fpc = 0.2) was injected to 
fit the Auger spectrum [52]. The source emissivity is assumed to be constant over time. Newly born 
pulsar population are assumed to have log and P normally distributed, and wind acceleration effi- 
ciency r] = 0.3. Simulation results were normalized at 10^^ eV with /g ~ 0.05. Notice the propagated 
spectrum and composition remain almost unchanged compared to the Auger-uniform case with 50% 
Proton, 30% CNO and 20% Fe injection (top panels in Fig 1 and Fig 2). 
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Figure 4. Cosmic ray flux measurements by KASCADE-Grandc [72], Auger [15] and TA [73] com- 
pared with pulsar model predictions. The total spectrum in solid black sums up cxtragalactic (dash) 
and Galactic (solid) components. The top, middle, and bottom panels correspond to the Auger- 
uniform, Auger, and TA cases respectively, as in Fig 1. Pulsar and propagation parameters: wind 
acceleration coefficient 77 = 0.3, Galactic magnetic field coherence length Ic = 20 pc, magnetic halo 
height H = 2 kpc. 
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Figure 5. Average logarithmic mass of cosmic ray derived from X^ax measurements from [75] 
with data from Tunka [76], Yakutsk [77, 78], CASA-BLANCA [79], HiRes/MIA [80], HiRes [81], 
KASCADE-Grande [72], Auger [15] and TA [73] for hadronic interaction model EPOS vl.99 [57] 
compare with simulation predictions (red lines) as in Fig. 4. Dashed lines indicate the energy range 
where pulsars contribute less than 80% to the total flux (see Fig. 4) and other sources also contribute. 
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